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ABSTRACT 

We present an analytical method for studying the changes of the orbital 
characteristics of binary systems with circular orbits due to a kick velocity 
imparted to the newborn neutron star during a supernova explosion (SN). 
Assuming a Maxwellian distribution of kick velocities we derive analytical 
expressions for the distribution functions of orbital separations and eccentricities 
immediately after the explosion, of orbital separations after circularization of 
the post-SN orbits, and of systemic velocities of binaries that remain bound 
after the explosion. These distributions of binary characteristics can be used 
to perform analytical population synthesis calculations of various types of 
binaries, the formation of which involves a supernova explosion. We study in 
detail the dependence of the derived distributions on the kick velocity and 
the pre-SN characteristics, we identify all the hmits imposed on the post-SN 
orbital characteristics, and we discuss their implications for the population 
of X-ray binaries and double neutron star systems. We show that large kick 
velocities do not necessarily result in large systemic velocities; for typical X-ray 
binary progenitors the maximum post-SN systemic velocity is comparable 
to the relative orbital velocity prior to the explosion. We also find that, 
unless accretion-induced collapse is a viable formation channel. X-ray binaries 
in globular clusters have most probably been formed by stellar dynamical 
interactions only, and not directly from primordial binaries. 



Subject headings: stars: binaries - stars: neutron - stars: supernovae - X-rays: 
binaries 
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1. INTRODUCTION 

Studies of the radio pulsar population (e.g., Gunn & Ostriker 1970; Helfand & 
Tademaru 1977; Harrison, Lyne & Anderson 1993; Lyne & Lorimer 1994) have shown that 
pulsars move in the Galaxy with very high space velocities, ranging from 20 to 2000 km s~^, 
and that their galactic distribution has a large scale height, of the order of Ikpc. The 
origin of these high velocities is often attributed to a kick velocity imparted to the neutron 
star at the time of the supernova explosion. Early studies by Dewey & Cordes (1987) and 
Bailes (1989) concluded that the mean magnitude of the kick velocity is of the order of 
100 — 200kms~^. In a more recent study, which takes into account new measurements of 
pulsar proper motions, a new electron density model, and a selection effect against fast 
pulsars, Lyne & Lorimer (1994) found the mean kick velocity to be ~ 450 ± 90kms~^. 
Additional observational evidence in support of a kick velocity imparted to neutron stars at 
birth are related to the existence of a high-velocity population of 0,B runaway stars (e.g.. 
Stone 1991), as well as to supernova remnant - pulsar associations, studies of which yield 
kick velocities up to 2000 km s"^ (Caraveo 1993; Frail, Goss, & Whiteoak 1994). 

In contrast to Dewey & Cordes (1987), Iben & Tutukov (1996) have recently concluded 
that the hypothesis of natal kicks imparted to neutron stars is unnecessary. They have 
found that the transverse velocity distribution of pulsars in the solar neighborhood, as well 
as that of O, B runaway stars, massive X-ray binaries, and double neutron stars, can be 
explained by the recoil velocity due to symmetric supernova explosions. However, they 
reach this conclusion by assuming (i) that all stars are members of binary systems and 
(ii) that neutron stars formed by massive single stars (formed only by mergers) or in wide 
binary systems rotate too slowly to become radio pulsars. Although, their results are 
marginally consistent (mean predicted velocities are ~ 100 — 150kms~^) with the old pulsar 
distance scale (Harrison et al. 1993), they are not consistent with the more recent results of 
Lyne & Lorimer (1994). 

Over the years, several theories have been put forward in an effort to explain the origin 
of kick velocities (e.g., Harrison & Tademaru 1975; Chugai 1984; Duncan & Thomson 1992; 
Herant, Benz & Colgate 1992; Janka & Miiller 1994; Burrows, Hayes, & Fryxell 1995; 
Burrows & Hayes 1996). Even a small asymmetry during the collapse of the core can give 
a kick to the remnant of the explosion. The asymmetry may be related either to neutrino 
emission or to mass ejection during the supernova, and may be caused by the magnetic field 
or rotation of the collapsing core, or by hydrodynamic instabilities, such as Raylcigh- Taylor 
or convective motions. In any case, the mechanism responsible for the kick velocity is still 
not well understood, and it appears that fully three-dimensional numerical simulations of 
the core collapse will be required in order to settle this issue. 
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Several authors have previously studied the effect of an asymmetric supernova explosion 
on binary parameters, focusing on various aspects of the problem. Early work by Flannery 
& van den Heuvel (1975), Mitalas (1976), Sutantyo (1978), and Hills (1983) addressed the 
problem of deriving expressions of post-SN orbital characteristics for a specific kick velocity 
for both circular and eccentric pre-SN orbits. They also derived survival probabilities for 
kick velocities of constant magnitude and random direction. The one-to-one link between 
pre-SN and post-SN parameters is broken when kick velocities are allowed to have a 
distribution over both magnitude and direction, in which case there exists a distribution of 
post-SN characteristics, even for pre-SN binaries with specific orbital parameters. Wijers, 
van Paradijs, & van den Heuvel (1992) were the first to address this problem. They derived 
an analytic expression for the distribution of post-SN orbital separations and eccentricities 
only, which however was also convolved with a distribution of pre-SN orbital separations. 
More recently Brandt & Podsiadlowski (1995) addressed the same problem using numerical 
methods (Monte Carlo simulations). The resulting distributions are again convolved with 
pre-SN period distributions, and, in this case, are calculated only for specific stellar masses, 
in an effort to compare them with observation. Because these distributions are calculated 
numerically, information about the allowed ranges of post-SN characteristics, the shape 
of multi-dimensional distributions, and their dependence on the pre-SN and kick-velocity 
characteristics is limited. Our purpose here is to derive analytical expressions of various 
post-SN characteristics for the realistic case of kick velocities with a distribution in both 
magnitude and direction. The derived distributions are general, and apply to any circular 
binary systems that experience asymmetric supernova explosions. 

The study presented in this paper has been motivated by our interest in performing 
population synthesis calculations for low-mass X-ray binaries. Monte Carlo techniques 
have been widely used in such calculations modeling various kinds of binary systems 
(e.g., Dewey & Cordes 1987; de Kool 1992; Romani 1992). Another method is based in 
creating a multi-dimensional grid of initial binary parameters and tracing the evolution of 
systems through a sequence of evolutionary stages for each set of initial parameters (Kolb 
1993; Iben, Tutukov, & Yungel'son 1995 and references therein). Both of these numerical 
methods have the same problem: although the goal is to calculate the characteristics of 
the final population, the sampling procedure is applied on the initial population, and 
therefore it is possible that the final population is under-sampled, even if the sampling 
of the primordial population appears adequate. Another problem with both methods is 
related to statistical accuracy: the initial sets of parameters cover a wide range, of which 
only a small part is populated by progenitors of interest, especially in the case of X-ray 
binaries, which have very small birth rates; therefore, with these methods it is necessary 
to study a very large number of primordial binaries in order for a statistically significant 
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number of systems to survive. Both of these problems are absent in population synthesis 
calculations performed analytically, where distributions of primordial binaries over orbital 
characteristics are transformed through a sequence of evolutionary stages, using Jacobian 
transformations. In this way, the regions in parameter space populated by the progenitors of 
interest are identified and the final population is calculated directly. This method has been 
formulated and applied to the study of cataclysmic binaries by Politano (1996) (see also 
Politano, Ritter, & Webbink 1989; Politano & Webbink 1990). Apart from the absence of 
the problems discussed above, the analytical method has the additional advantage that the 
shape of final distributions is calculated exactly, revealing fine details and subtle features, 
such as sharp peaks, infinities, or definite limits imposed on the final parameters. Also, the 
various dependences of these parameters and their distributions on the initial parameters 
can be identified and studied in detail. 

In order to use the analytical method in population synthesis of neutron-star binaries, 
it is necessary to develop an analytical tool for the modeling of asymmetric supernova 
explosions. Our purpose in this paper is to present such a method based on Jacobian 
transformations for computing analytically the probability distributions of several orbital 
characteristics of post-SN binaries. These distributions include (§2) the orbital separations 
and eccentricities immediately after the supernova explosion, the circularized orbital 
separations, and (§ 3) the systemic velocities. They are derived for circular pre-SN orbits 
and for kick velocities that are randomly distributed not only in direction but also in 
magnitude (Maxwellian distribution). The derived expressions can be used in synthesis 
calculations of any kind of binaries that experience supernova explosions during their 
evolution. The analytical character of the derivation enables us to perform detailed 
parameter studies and identify those characteristics of the kick velocities or the prc-SN 
binaries that govern the behavior of the post-SN distribution functions. Expressions for 
the limiting cases of very large or very small kick velocities relative to the prc-SN orbital 
velocities are also derived. We identify the limits imposed on the post-SN parameters and 
discuss their physical interpretation. In addition, we calculate survival probabilities (§ 4) as 
functions of the pre-supernova (pre-SN) orbital characteristics and the mean kick velocity. 
Finally, we examine several implications of our results (§5) for the progenitors of high- and 
low-mass X-ray binaries, double neutron stars, and their populations in globular clusters. A 
list of the symbols used throughout the paper is given in Appendix A. The study of a few 
special cases is included in Appendices B and C. 
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2. POST-SUPERNOVA ORBITS 

We assume that the binary orbits prior to the supernova explosion are circular, and 
that the kick velocities follow a Maxwellian distribution. The first of these two assumptions 
is unlikely to be violated in the case of systems that have emerged from a common-envelope 
phase, or have experienced (semi-)conservative mass transfer. Orbital eccentricities may 
be important for binaries with components that have not interacted prior to the supernova 
explosion, although it is still possible that circularization has occurred during their 
main-sequence evolution (see Portegies-Zwart & Verbunt 1996). The second assumption, 
that of a Maxwellian distribution of kick velocities, we adopt in the absence of an adequate 
theoretical understanding of their origin. It is conceivable that the direction of the kick 
is affected by the kinematical or rotational properties of the collapsing core, but any 
correlation between the kick direction and the orbital rotational axis or the orbital velocity 
has yet to be established. In principle, a method like that described here can be used with 
any distribution of kick velocities, although the derivational details will be different. 

In most of our calculations the interaction between the expanding supernova shell and 
the companion to the exploding star has been ignored. According to Fryxell & Arnett 
(1981), this interaction is generally weak, especially in the case of low-mass X-ray binaries 
(LMXBs) and double neutron stars, where the solid angle intercepted by the companion is 
very small unless the orbital separation is also small (see Appendix A and Romani 1992). 
This effect may be more important for high-mass X-ray binaries (HMXBs), although their 
orbits are much wider. Nevertheless, for completeness we have repeated the calculation of 
the probability distribution of circularized orbital separations, including also the effect of 
the impulse velocity (see Appendix A). 

2.1. Non- Circularized Orbits 

In this section we derive the distribution of post-SN binary systems over orbital 
separations and eccentricities immediately after the supernova explosion, by transforming 
the distribution of kick velocities imparted to the neutron star into a distribution over other 
binary parameters of interest. 

We adopt a reference frame centered on the exploding star mass. Mi, just prior to the 
supernova. The companion mass, M2, is chosen to be at rest and the exploding star to move 
in a circular orbit with separation Ai. The x-axis lies along the line connecting the centers 
of mass of the two stars, pointing from M2 to Mi. The y-axis lies parallel to the direction 
of the pre-SN orbital velocity K of Mi relative to its companion. The z-axis completes a 
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right-handed orthogonal system (Figure 1). 

Two of the parameters characterizing the post-SN binaries are the orbital separation, 
Af, and the eccentricity, e. We use the energy and angular momentum equations for 
eccentric orbits to relate these two parameters to the three components of the kick velocity. 
The supernova explosion, mass loss, and kick are assumed to be instantaneous. In general, 
for two stars with masses Ma and Aff, in an orbit with orbital separation A and eccentricity 
e, their relative velocity at a distance r is given by: 

V = G{Ma + M,) (1) 

and the specific angular momentum of the system is: 

\fxV\^^G{Ma + Mb)A{l-e^). (2) 

We can also apply the above two equations in the case of the post-SN binary, for which 
V — {Vkx: Vky + Vkz) ^ud v — (Aj, 0, 0). The orbital separation, Af, and the eccentricity, 
e, of the post-SN orbit are thus related to the components of the kick velocity, Vkx, Vky, and 
Vkz by the expressions: 



G{Mns + M2) 



2G{Mms + M2) 



Ai 



-1 



(3) 



, 2 ^ + Vky + V? + '^VkyVMj 

^ GiMMs + M2)Af ' ^ ' 

where M^s is the gravitational mass of the neutron star and Vk is the magnitude of the 
kick velocity. 

The third independent parameter describing the post-SN state of the binary is the 
orientation of the eccentric orbit relative to the pre-SN orbital plane. The plane of the 
binary orbit is altered due to the z-component of the kick velocity. Since the explosion is 
assumed to be instantaneous, the position of the two stars just before and just after the 
supernova remains unchanged. In our reference frame (Figure 1) the two stars lie along the 
X-axis, and therefore the intersection of the two orbital planes must coincide with the x-axis. 
The angle, 9, between the pre- and post-SN orbital planes, is equal to the one between the 

— * 

relative velocity just before the explosion, Vr — (0, K, 0), and the projection of the relative 
velocity just after the explosion onto the y — z plane, Vy^ — (0, Vky + Vy, Vkz)- Hence: 

cose = = + (5) 
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For convenience, we rewrite equations (3), (4), (5) in dimensionless form using the 
following definitions: 



a 



P = 



A' 

Mns + M2 
M1 + M2 ' 



(6) 
(7) 
(8) 



where j can be x, y, or z. The dimensionless post-SN orbital separation a, the eccentricity 
e, and the angle 9 between pre- and post-SN orbital planes can be expressed as functions of 
Vkx, Vky, Vkz, and /3: 



a — 



/3 



1 - e 
cos^ = 



Hz + (t'/c. + 1)^] 



1/2- 



(9) 
(10) 
(11) 



To obtain the distribution of binaries over post-SN characteristics we use the Jacobian 
transformation of the kick velocity distribution. The distribution functions for each of the 
three components of the kick velocity are assumed to be Gaussian: 



(12) 



where ^ = a/Vr, and a is the velocity dispersion of each of the one-dimensional Gaussian 
distributions. Using equations (9), (10), and (11) we obtain: 



g{a, e, cos^) 



P 



3/2 



2 e 



\2ne 
X exp 

X exp 



2^lV2 



2e 



[a{l - e^)] 



' ^/5^ + l 



a — 



1 + eJ Vl-e 



— a 



-1/2 



a 



2^^l/2 



■ cos^ 



(1 - cos^ e) 



-1/2 



(13) 



In general, we are not interested in the orientation of the post-SN orbital plane. 
Therefore, we integrate 51(0;, e, cos ^) over all orientations, and obtain for the distribution 
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over orbital separations and eccentricities: 

3/2 



G{a, e) = 




27ie 



1 



1 2a- 1 ^ 
7^/5 + 1 



a 



2^ 



a 



1 + e 



1-e 



-1/2 



(14) 



where 



(/3 « (1 - e^)) 



1/2 



and Jo is the modified Bessel function of zeroth order. The above expression has two 
singularities, at a = 1/(1 ± e), which correspond to the special case of the velocity of the 
newborn neutron star being restricted in the y-z plane, V^x — 0. We can see this by using 
equations (9), (10) to obtain: 



P 



a 



a — 



1 + e 



a 



(15) 



In the singular cases, the distribution of kick velocities becomes two-dimensional, and there 
arc only two independent variables describing the post-SN state: cos^ and e (or a). The 
corresponding distributions arc derived in Appendix B. 

Scrutiny of equation (14) indicates that post-SN binaries populate only a restricted 
area of the a — e plane, independent of the orbital characteristics of the pre-SN systems. 
Acceptable values for a span a range from 1/(1 -|- e) to 1/(1 — e), limits which were first 
identified by Flannery & van den Heuvel (1975). Since the post-SN orbit must include the 
position of the two stars just prior to the explosion, the post-SN orbital separation, Af, 
cannot be smaller than half of the pre-SN separation, A^. 

In addition to remaining bound, the two stars in the post-SN binary must avoid 
physical coUision, which would probably lead to a merger. The closest distance between the 
two stars must at a minimum exceed the sum of their radii: 



or 



Af (1-e) > Rns + R2^R2, 



a (1 — e) > — ^ = c. 

^ ^ Ai 



(16) 



(17) 



This condition sets a lower limit on a, o; > c/(l — e), or an upper limit on e, e < 1 — c/a. 
The complete set of limiting curves on the a — e parameter space, for a range of different 
values of c, is shown in Figure 2. Nevertheless, we will set c = for simplicity in the 
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foUowing discussion, returning at the very end of this section to comment on the effect of a 
non-zero value of c. 

A two-dimensional distribution over a and e (eq. [14]) is shown in Figure 3 for the 
specific choice of /3 = 0.6 and ^ = 1. The behavior of the distribution is dominated by the 
square root term that appears in equation (14). This term becomes equal to zero along the 
q;(1 ± e) = 1 curves. Variation of the values of (5 and ^ affects only the normalization of the 
distribution, and not its qualitative shape. 

The distribution of post-SN systems over eccentricity, J{e), can be found by integrating 
G(q;, e) over a. Sample distributions of post-SN systems over e are plotted in Figures 4a 
and 4b for different values of j3 and ^. For /3 > 0.5 (Figure 4a), less than half the total 
mass of the pre-SN system is lost, and the binary would remain bound in the case of a 
symmetric explosion (no kick imparted to the neutron star). In the limit that ,^ ^ 0, the 
distribution over e sharply peaks at e — > {1 — (3)/(3, which is the eccentricity of the post-SN 
orbit if the explosion were symmetric (see, for example Verbunt 1993). As ^ increases, and 
the kick velocity becomes comparable to the relative orbital velocity of the stars prior to 
the supernova explosion, the distribution becomes broader, and then declines uniformly for 
^ 1. When /3 < 0.5 (Figure 4b), binaries would be disrupted in the absence of any kick 
velocity. In the limit that ^ — > 0, few systems survive generally with very high eccentricities 
(e 1). As ^ increases, JT{e) grows until ^ ~ 1 — ^/J3, then declines. For ^ S> 1, i7(e) 
converges to the same asymptotic form regardless of whether beta > 0.5 or not: 

where 

/ 2e 

P = 



. 1 + e 

and K{p) is the complete elliptic integral with the following series representation: 

1 2 



if(p) = | l + (i)V + (^)V + ...+ 



(2n - 1) 



where n is a positive integer. 

To derive the distribution over orbital separations, Q{a), we integrate over eccentricities, 
e. Plots of the distribution over a are given in Figures 5a and 5b for a set of different 
P and ^. When (3 > 0.5 and ^ — > (Figure 5a), conditions are similar to those in a 
symmetric supernova. The distribution is narrow and peaks at a — > /3/ (2/3 — 1), which is the 
post-SN orbital separation in the absence of any kicks (see Verbunt 1993). As ^ increases 



the distribution broadens and peaks at orbital separations smaller than that before the 
explosion. For /? < 0.5 (Figure 5b) and small kick velocities, the probability of disruption 
is very high and the binaries that survive have very large orbital separations. For higher 
values of ^, more systems are able to reduce their energy and remain bound with smaller 
orbital separations. When ^ ^ 1, Q{a) converges to the asymptotic form: 



where a > 1/2. We note that for ^ 1, both Q{a) and i7(e) assume forms which, apart 
from their normalization do not depend on either ^ nor on /3. 

In the preceding discussion, we have not accounted for the possibility that the two 
stars may collide after the supernova explosion. As a consequence of this last constraint the 
parameter space in a and e is further restricted, and the integrated distributions of orbital 
separations and eccentricities are altered (see Figures 6a and 6b). It is evident that the 
survival probabilities decrease dramatically if c is not very small (c ^ 0.01), especially for 
large eccentricities and for orbital separations larger than that of the pre-SN binary. 



Tidal interaction between the binary members leads to circularization of the post-SN 
orbit, on a time scale that depends on the characteristics both of the eccentric orbit and of 
the companion to the neutron star. Setting aside the question of the relevant time scales 
we can calculate the distribution of post-SN systems over orbital separation, after 
circularization has been achieved. 

During the circularization process, orbital energy, E, is dissipated while orbital angular 
momentum, J, is conserved. We define the dimensionless quantities: 




(19) 



2.2. 



Circuleirized Orbits 




(20) 



E 



1 



(21) 



e 



Eo 



a 



where Jl = G AiM%sM^/{MNs + M2) and Eo = -G MnsM2/ {2Ai). From conservation of 
orbital angular momentum we find for the orbital separation, A^ of the circularized orbit: 
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Using equations (20), (21), and (22) we transform G{a,e) into the distribution of post-SN 
circularized systems over etc and e: 



H{ac, e) — n 



/3 



3/2 



exp 



2/3 + r 

2^ ^ 



exp i^^e 
V2 - die - e ' 



(23) 



where (2 — etc — e) > 0. 



The circularized post-SN systems are characterized by only one parameter, the orbital 
separation. In order to find their distribution over ac we need to integrate H{ac, e) over the 
dimensionless orbital energy e. The limits of integration are found by considering all the 
constraints that viable post-SN systems must satisfy. 

An upper limit to e is set by the geometrical constraint, that the post-SN eccentric 
orbit must include the position of the stars prior to the supernova explosion (see eq. [23]), 

e < 2 -etc. (24) 

The second constraint is that the post-SN system must be bound, and hence its orbital 
energy E must be negative. Since Eo has been defined to be negative, we obtain a lower 
limit for e: 

e > 0. (25) 

An additional lower limit is set by the need to avoid a physical collision. This condition is 
expressed as a lower limit on a, i.e., a > c/{l — e), such that the periastron distance in the 
eccentric orbit exceeds the radius of the companion to the neutron star. Using equations 
(20), (21), and (25), we can rewrite this condition as: 

e > ^. 26 



By checking for consistency {emax > ^min) among the above limits, we find that from 
all possible values of etc only a small range is acceptable for post-SN circularized orbits: 

< ac < 2. (27) 
1 + c ^ ^ 



The post-SN systems may be divided into groups depending on the value of etc- 
Systems with etc > 2 become unbound (e becomes negative). Systems with etc < 2c/ (1 -|- c) 
are bound, but all lead to a merger of the two stars. For 2c < etc < 2, systems are bound 
and they all avoid collision. In this case, e^j„ = and e^ax — 2 — ac- Finally, systems with 
2c/(l + c) < ac < 2c are also bound, but a fraction of them merge. For this reason the 
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range of acceptable energies is further constrained: e^m = (2c — ccc) /c^ > 0. Clearly, unlike 
the case of a symmetric explosion, where the STipcrnova always rcsTiIts in an expansion 
of the orbit (Verbunt 1993), the separation of the circularized orbit after an asymmetric 
explosion may become smaller than the pre-SN separation. However, although the lower 
limit on a is extended to values smaller than unity due to the kick velocity, the upper limit 
to the post-SN circularized orbital separation remains twice the prc-SN separation. 

To obtain the distribution of post-SN binaries over etc we integrate over e. This 
integration can be performed analytically, yielding: 

where: 

erf(xo) 



^ / e dx, 

VTT JO 

r 2c- 

V2-«c 



2c 



< ac<'2,c 



= V2^^, 



2c < ttc < 2. 



The behavior of H{ac) is dictated by the values of the two parameters (3 and ^. Using 
the asymptotic forms of the error function and the modified Bessel function in the two 
limits that the r.m.s. kick velocity is much larger or much smaller than the relative orbital 
velocity in the pre-SN orbit we obtain: 



lim 7Y( etc) 



1 fp- 



3\ 1/4 



2\/27r V«c 



C ^exp 



3/2 



(29) 
(30) 



The distribution over ac for different values of f3 and C, is plotted in Figure 7. The behavior 
of Ti.{ac) is analogous to that of G{a, e). In the limit of kicks much smaller than the relative 
velocity of the stars in the pre-SN orbit, conditions approximate the case of a symmetric 
explosion, and the distribution peaks at those values of ac that are consistent with such an 
explosion: etc = 1//? for /? > 0.5 (Verbunt 1993), and etc = 2 for /3 < 0.5. As the average 
kick velocity becomes comparable to K, smaller orbital separations become more abundant, 
and for even larger kicks, a shrinkage of the orbit relative to the pre-SN state is favored. 
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3. SYSTEMIC VELOCITIES 

During the supernova explosion, the post-SN system as a whole receives a velocity 
relative to the center of mass of the pre-SN binary. We derive the probability distribution of 
the systemic velocities by performing a sequence of Jacobian transformations of the initial 
distribution of kick velocities. 

We choose to work in a reference frame, in which the three axes have the same 
orientation as the one shown in Figure 1, but which is centered on the center of mass of the 
system prior to the supernova explosion. In this frame the vector velocities of the two stars 
are: 

Mo 

where Vr is the magnitude of the relative velocity of the two stars. After the supernova 

— * — * 

explosion, V2 remains the same and Vi becomes: 

Vns = {Vk,, Vky + Vr, V,,). (33) 

Ml + M2 

Hence, the systemic velocity is: 

^ ^ MnsVns + M2V2 
Mms + M2 

= ^ (MNsVk^, MusVky - ^^\~^^fJ^^ Vr, MNsVk,] . (34) 

MNS + M2 V M1 + M2 J ^ ' 

We define the dimensionless systemic velocity Vgys = V^ys/K- Using equations (9), (10), 
and (11), we obtain: 



2a 



V = + 1^3 cos 9 a{l - e ) , (35) 



1/2 



where: 



Ki = 



(Mi + M2)2' 

K2 = 



Mis 



(M^5 + M2)(Mi + M2)' 
K3 = IsJkx K2. (36) 
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We derived above an expression (eq. [13]) describing the distribution of post-SN systems 
over orbital separation, a, eccentricity, e, and orientation of the orbital plane, cos^. Using 
equation (35) we can eliminate cos^ and transform g{a,e,cos9) into a distribution over a, 
e, and the magnitude of the systemic velocity, Vsys' 



S\Cki 6, Vsys ) 




a — 



1 + e 



1-e 



— a 



-1/2 







exp 





sys) 



21 



1^1 + H2- 
1/2 



2a -I 



a 



^sys 



(37) 



The above expression is valid in the general case of v^^ 7^ and v^z 7^ 0. The special case 
of Vkx — has already been discussed and corresponds to the pre-SN orbital separation 
becoming either the periastron or the apastron distance in the post-SN eccentric orbit (see 
eq. [15]). In the special case of v^z = 0, the plane of the orbit remains unaffected by the 
explosion, since the kick velocity is restricted in the x — y plane, which is the orbital plane 
prior to the explosion (see eq. [11], [35] and Figure 1). The derivation of the probability 
density for these special cases is described in Appendix B. 

It is also interesting to study how the systemic velocity imparted to the binary during 
the supernova explosion correlates with the orbital separation after the circularization. We 
transform s{a, e, Vgys) to a distribution over {ac, e, Vsys) (see eq. [20], [21], [22]): 



f{ac,e,Vsys) 




XV 



sys 



X 



(2 — ttc — e) ^^"^ exp 

(^Kl + 2k,; 




^ K2e + Viys 



K2€ — V 



sys 



-1/2 



(38) 



and we need to integrate over orbital energies, e. The above distribution has three poles in 
e; for clarity we rewrite it as: 



p Y^' 2 



sys) 



K,2 



exp 



2e 



2^ V 



\ (2;^^^^ -(2/3 + 1) 
2 



«3 



XV sys exp 

X(6-A0-^/^ (A2-6)-V^ (A3-6)-V^ 



(39) 
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where 



A 



1 




A3 



2-ac- 



All three poles are numerically integrable except in the special case that: 



A2 = A3 ^ V, 



sys 



(40) 



In this case, the two-dimensional distribution FioLc^Vgy^ becomes infinite along the line 
defined by equation (40). However, it is still integrable over and Vsys^ so that the final 
integral is finite. A sample distribution Fiptc^Vsys) for a specific choice of Mi, M2, and 
and for c = is shown in Figure 8. The spikes correspond to the pole (eq. [40]). The limits 
on ccc are in agreement with equation (27), while the limits on Vgyg are discussed in the next 
section. 

We can obtain the distribution of systemic velocities only, J-'{vsys), by further 
integrating over etc- Sample distributions normalized to the total survival fraction for 
different values of ^ are given in Figure 9. It is evident that the limits imposed on Vgys 
are independent of the characteristics of the kick velocity distribution. As the magnitude 
of the kick velocity increases the systemic velocity remains restricted to a certain range 
of values specified by the stellar masses. Within this range, the distribution function of 
Vsys shifts towards larger velocities as ^ increases, reaching an asymptotic distribution for 
^ ;^ 3, independent of ^ (sec Figure 9). As in the case of the one-dimensional distributions 
over eccentricities and orbital separations (see eq. [18] and [19]), this behavior is due to the 
fact that for ^ ^ 1, the exponential terms in equation (38) approach unity, and only the 
normalization constant depends on ^. In this limit the normalized distribution of systemic 
velocities depends only on the stellar masses involved. 



Scrutiny of equation (38) shows that there exists an upper and a lower limit to the 
values of the systemic velocities, since the expression shown is real only if 



3.1. Limits on the systemic velocity 



KgOic - (ki + 2^2 - «2e - Vgys)"^ > 0) 



(41) 



or 



{ki + 2K2 - K,2e - K3y/a^)^^^ < Vsys < («1 + 2^2 " + K3"\/c^)^* 



(42) 
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In addition, we have already derived limits on e and etc (eq. [24] - [27]). By taking these into 
account, we can find the absolute lower and upper limits on Vgys'- 

Vsys < \f^\ + V^2k^ 

Vsys > |a/^ - \/«cK2| > x/kT - V2k^, if — > 2. (43) 

The last inequality is true if Mi > 2Mjv5 = 2.8 Mq (eq. [36]), assuming a neutron star 
gravitational mass equal to 1.4 Mq. This condition is satisfied for all the progenitors of 
LMXBs forming via the He-star and direct supernova mechanisms (Kalogera & Webbink 
1996; Kalogera 1996), and for most of the HMXB progenitors (Portegies-Zwart & Verbunt 
1996). Therefore, the maximum and minimum systemic velocities are: 

Ml + M2 {Mns + M2)V2(Mi + M2)V2 ^ > 

v^P = ^^^\^ - jjj .i^f/o){! TTvTT?' foi^ Ml > 2Mns (45) 

M1 + M2 (M^5 + M2)V2(Mi + M2)V2' V y 

The above upper limit on Vgys is in agreement with the one found by Brandt & Podsiadlowski 
(1995), while the lower limit (eq. [45]) is stricter than theirs. Inspection of equation (44) 
shows that a maximum value of v'^f^ exists in the limit that (i) M2 is equal to zero and 
(ii) Ml is equal to the minimum possible mass of a neutron star progenitor, ~ 2.2 Mq 
(e.g., Habets 1985). In this limit, the maximum of f™"^ is ~ 2, and hence the systemic 
velocity of a bound post-SN binary can never exceed twice the value of the prc-SN relative 
orbital velocity, regardless of the magnitude of the kick velocity and the masses involved. 
Therefore, it becomes clear that high kick velocities do not necessarily result in high 
systemic velocities, as well. 

Both upper and lower limits on systemic velocities can be understood physically. The 
maximum systemic velocity is acquired by that binary for which the neutron star receives 
a kick velocity oriented opposite to the pre-SN orbital velocity with a magnitude, such 
that its post-SN kinetic energy is just below its binding energy. The minimum systemic 
velocity is acquired by that binary in which the neutron star receives a kick velocity with 
the smallest possible magnitude needed to avoid disrupting the system due to mass loss. It 
is important to re-emphasize that neither the upper nor the lower limits depend on the kick 
velocity distribution. 



4. SURVIVAL PROBABILITIES 



The total survival probability of a binary system with specific initial orbital 
characteristics is of interest to studies of the statistical properties of an entire population 
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of binaries. We can obtain this survival probability by integrating over the distribution 
of circularized dimensionless orbital separations, T-C{ac)- This integration can only be 
performed numerically, but it is straightforward since the function has no poles and the 
limits of integration are well-defined (eq. [27]). 

Clearly, kick velocities to neutron stars may bind systems that would otherwise be 
disrupted, or disrupt those that would have remained bound. If the average kick velocity is 
large compared to the initial relative orbital motion of the binary components then survival 
depends on the small probability that the kick velocity is itself small and directed opposite 
to the original motion of the collapsing component. If the ratio of kick velocity to initial 
relative orbital velocity is small, then the survival rate of systems that would otherwise be 
disrupted falls very rapidly as this ratio decreases. Asymptotically, we have, respectively 



The survival fractions for a wide range of values of ^, and for two different values of (3, 
are shown of Figure 10. In this illustration, ^ is varied by keeping the r.m.s. kick velocity 
constant, while allowing the pre-SN orbital separation to vary. Among those systems 
that would remain bound if collapse were symmetric (Mi = 3.8 Mq and M2 — 1.0 Mq), 
kick velocities will tend to unbind widely separated binaries, for which the relative orbital 
motion falls below the kick velocity. Among systems that suffer so much mass loss in a 
supernova that they would otherwise be disrupted (Mi = 8.6 Mq and M2 = 1.0 Mq), kick 
velocities will favor the survival of binary systems in just that range of separations where 
the relative orbital velocity is comparable to the mean kick velocity. 



The expressions derived here provide a tool necessary in analytical population syntheses 
of neutron star binaries. The additional step needed in such syntheses is to convolve 
the distribution over post-SN parameters with the distribution of pre-SN binaries over 
masses and orbital separations (see also Wijers et al. 1992). This link depends on the 
type of final systems and the specifics of their formation mechanism. In addition, the 
distribution functions of systemic velocities and their correlation with orbital separations 
and eccentricities (or circularized orbital separations) can be used in studying the motion 




(46) 



5. 



CONCLUSION 
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of neutron star binaries in the Galactic potential, and in modeling their spatial distribution 
in the Galaxy. 

The results of the study presented in this paper have a number of important 
implications concerning the population of neutron star binaries: 

There exists a correlation between orbital separations and eccentricities, which is 
independent of the characteristics of the binary or the magnitude of the kick velocity. 
For post-SN orbits much wider than the pre-SN orbit, the total energy of the binary 
significantly increases, and the system remains bound only in a highly eccentric orbit. On 
the other hand, the eccentricity may be low (e ^ 0.4) only if the post-SN orbital separation 
is comparable to that before the explosion. The discovery of a double neutron star system 
of modest eccentricity could therefore be used to infer the size of the orbit of its progenitor, 
provided that the gravitational radiation decay time scale for the orbit were long enough 
for such losses to be negligible. 

The ratio of the post-SN systemic velocity, Vgys, to the pre-SN relative orbital velocity, 
Vr, is restricted in a relatively narrow range of values. Both lower and upper limits depend 
only on the stellar masses involved. For the ranges of progenitor masses relevant to HMXBs, 
LMXBs, and double neutron star binaries we find V^™"^ ^ 1.5 Since LMXB progenitors 
are more tightly bound than those of HMXBs, and hence have higher relative orbital 
velocities than HMXB progenitors, the systemic velocities of LMXBs are expected to be 
higher than those of HMXBs. It is also clear that measurements of systemic velocities 
of neutron star binaries do not necessarily reveal information about the kick velocities 
imparted to neutron stars in individual systems. Instead, they can be used to infer typical 
relative orbital velocities prior to the supernova explosion. 

Although the allowed range of systemic velocities is independent of the kick velocity, 
the probability distribution within this range does depend on the r.m.s. magnitude of 
kick velocities relative to the pre-SN orbital velocities. In the limit of very small kicks 
the distribution sharply peaks at values close to the lower end of the range. As the r.m.s. 
of the kick velocities increases the distribution becomes broader and its peak shifts to 
higher velocities. For kicks much higher than the pre-SN relative orbital velocities, the 
shape of the distribution remains unaffected, and further increases of the r.m.s. kick 
velocity only decrease the binary survival rate, without altering the velocity distribution of 
bound post-SN systems. Measurements of the systemic velocities of neutron star binaries 
will possibly prove quite significant in distinguishing between symmetric and asymmetric 
explosions with high kicks imparted to the neutron stars. 

The incidence of X-ray binaries in globular clusters relates to their smallest possible 
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systemic velocity and to how this velocity compares with the escape velocity from the 
cluster. For LMXBs formed via the explosion of the He-star remnant of a common envelope 
phase, typical parameters for the progenitors yield V^l'^ ~ 100 km s^^ (Kalogera & Webbink 
1996). The direct-SN channel (Kalogera 1996; Kalogera 1996) is fed by binaries with orbits 
which are much wider, but still small enough to avoid disruption by dynamical interactions. 
Typical parameters in this case yield VJ^^"' ~ 20 km s~^. Estimates of the escape velocities 
from the cores of globular clusters that contain LMXBs range from 30 km s~^ to 60 km s~^ 
(for NGC 1851, 6440, 6441, 6624, M15, and Lil 1); more loosely bound clusters such as Ter 
1 and 2 have central escape velocities of the order of lOkms"^ (Webbink 1985; van Paradijs 

1995) . It is therefore clear that post-SN binaries formed in globular clusters from primordial 
binaries via the He-SN channel have a very small chance of remaining in the clusters and 
becoming X-ray binaries. LMXBs formed via the direct-SN channel, on the other hand, 
will remain in the clusters, but their formation rate is too low to account for a significant 
fraction of the LMXB population in globular clusters. Barring accretion-induced collapse 
as an alternative formation channel, it therefore appears that low-mass X-ray binaries 
observed in globular clusters must have formed through stellar exchanges and captures, 
rather than directly from primordial binaries. 

We have already applied the analytical method presented here to study low-mass 
X-ray binaries formed via different evolutionary channels. The results of these population 
synthesis calculations will be presented elsewhere (Kalogera & Webbink 1996; Kalogera 

1996) . 

I am grateful to Ron Webbink for initiating this project, for many helpful discussions 
during the course of this study, and for carefully reading the manuscript. I would also like 
to thank Dimitrios Psaltis and Fred Lamb for valuable comments and discussions. This 
work was supported by National Science Foundation under grant AST92- 18074. 

Appendices 

A. Notation 

The symbols of the most important physical parameters are: 
Ml. mass of the exploding star, 

M2: mass of the companion to the exploding star and to the neutron star, 
Mjvs: mass of the neutron star. 
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i?2: radius of the companion to the neutron star, 

ol: ratio of the post-SN orbital separation to the pre-SN separation, 

e: eccentricity of the post-SN orbit, 

Q: angle between the pre-SN and post-SN orbital planes, 

/3: ratio of the total mass after the explosion to that before, 

^ : ratio of the standard deviation of the kick velocity distribution to the relative orbital 
velocity before the explosion, 

c: ratio of the radius of the companion star to the pre-SN orbital separation, 

ac- ratio of the circularized orbital separation to that before the explosion, 

e: ratio of the post-SN binding energy to that of a system consisting of the neutron star 
and the companion star in an orbit with the pre-SN orbital separation, 

Vsys- ratio of the systemic post-SN velocity to the pre-SN relative orbital velocity. 

The symbols of the various distribution functions over dimensionless orbital parameters, 
and respective equations, in which the corresponding expressions are given, are: 

^'(q;, e, cos^): distribution function of post-SN orbital separations, eccentricities, and cosine 
of the angles between pre- and post-SN orbital planes; equation (13), 

G{a^e): distribution function of post-SN orbital separations and eccentricities; equation 
(14), 

J{e): distribution function of post-SN eccentricities for ^ ^ 1; equation (18), 

Qipi): distribution function of post-SN orbital separations for ^ ^ 1; equation (19), 

HioLcie): distribution function of circularized orbital separations and post-SN binding 
energies; equation (23), 

T-C{ac): distribution function of circularized orbital separations; equation (28), 

s{a, e, Vsys)'- distribution function of post-SN orbital separations, eccentricities, and systemic 
velocities; equation (37), 

f{otc-i^-iVsys)'- distribution function of circularized orbital separations, post-SN binding 
energies, and systemic velocities; equation (38). 
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B. Effect of an impulse velocity 



According to the numerical calculations performed by Fryxell & Arnett (1981), 
approximately half of the momentum carried by the ejecta intersecting the companion is 
transferred to it. If Esn is the kinetic energy of the ejecta and Vimp is the velocity imparted 
to the companion of mass M2, then: 



MoV, 



1 



imp 



(Ml - M^s) 



2E. 



SN 



(Bl) 



where R2 is the radius of the companion. 



We assume that the impulse velocity, Vimp-, is given to the companion in a direction 
along the line connecting the two stars, pointing away from the neutron star. In the 
reference frame shown in Figure 1, the velocity of the neutron star relative to that of 

— * 

the companion immediately after the explosion is then: V — {Vkx + Vimp, Vky + K, Vkz)- 
Following the same procedure as that described in §2, we calculate the distribution of 
circularized orbital separations: 



T~iimp{oic) — 



(27re)3/^ 



exp 



2e 



^1, 



(B2) 



where 
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Ye 



erf 



— erf 



z± = V, 



imp ^ 13 [2 - ac - 

Vimp ± Y^/3(2 - Qic), 



2c — a. 



; 2c 

if < etc < 2c 



1 + c 
if 2c < etc < 2 



(B3) 



and Viinp = Vimp/Vr, Vr being the relative orbital velocity of the pre-SN binary. 

Integrating the above distribution over ac, we find the survival probability when the 
impulse velocity is taken into account. In Figure 11 we plot the survival probability as a 
function of the pre-SN orbital separation, Ai, with and without the effect of the impulse 
velocity, for the specific choice of Mi — 4Mq, M2 = 1 Mq (typical of an LMXB progenitor). 



and < >V2 



450kms~^. It is clear that the survival fraction decreases due to the 
impulse being imparted to the companion only when the orbital separation is very small, 
Ai ^SRq. This separation is smaller than typical values of the pre-SN orbital separations 
of progenitors of X-ray binaries or double neutron star systems. 
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C. Distribution over eccentricities and orbital separations for special cases 

C.l. Vk. = 

Following the same procedure as in the general case described in §2.1, we transform 
the two dimensional distribution of kicks Pyz{vky,Vkz) into g''(e,cos^), and then integrate 
over cos^. The eccentricity is given by: 



- 1 



(CI) 



where the plus sign corresponds to a{l — e) = 1, the minus sign to a(l + e) = 1, and cos^^ 
is still given by equation (11). The form of the integral of g'{e, cos 6) over cos^ is the same 
as in the general case, so we obtain: 



J'(e) = j;(e) + j:(e). 



(C2) 



where: 



J±ie) 



1 

2e 



exp 



;i±e)] 



-^(l + /3(l±e)) 



Io{Ze), 



1/2 



(C3) 



The plus sign again corresponds to being the periastron distance, and the minus sign to 
Ai being the apastron distance, in the post-SN orbit. The corresponding distribution over 
a is: 



2e oi 



2 ^2 



2e 



1 / 2a - 1 
1 + /3' 



-^o(^a) ) 



_ [/3(2a-l))/a]V2 



(C4) 



C.2. = 

In this case it is cos^ = ±1 (eq. [11]), that is, the post-SN orbital plane is either the 
same as the pre-SN one (^ = 0), or it has been rotated by an angle 6 = n, and the stars orbit 
in a retrograde sense after the explosion. We transform the two-dimensional distribution 
of kicks Pxy{vkx,Vky) into a distribution of dimensionless orbital separations and systemic 
velocities, s'{a,Vsys)- Using equations (9), (10), and (35) we obtain the expressions relating 
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the four variables: 



2 ^ 2a-l (3 _ 



ky 
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(C5) 
(C6) 



where the plus sign corresponds to cos^ = —1 and the minus sign to cos^ = 1. We calculate 
the necessary Jacobian and find: 



S'{a, Vsys) = S'_^_{a, Vsys) + S'_{a, Vgys), 



(C7) 



where: 
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Figure Captions 

Fig. 1 — Geometry of the binary system and reference frame adopted in the calculations. 
The orbital plane of the pre-SN binary coincides with the plane of the page, that is the x-y 
plane. 

Fig. 2 — Limits on the parameter space a — e of post-SN systems for a range of values of c, 
the ratio of the radius of the companion to the orbital separation prior to the explosion. For 
point stars {c — 0), the allowed parameter space is restricted between the two thick lines, 
which correspond to limits due to the geometrical constraint 1/(1 — e) > a > 1/(1 + e). The 
thin lines, which correspond to the constraint barring a physical collision, a > c/{l — e), 
impose a second, more stringent, lower limit to a. 

Fig. 3 — Distribution of post-SN systems over dimensionless orbital separations a and 
eccentricities e, for P — 0.6 and ^ = 1.0. 

Fig. 4 — Distribution of post-SN systems over eccentricities e for systems that (a) would 
remain bound (/3 = 0.6), or (b) be disrupted (/? = 0.4), in the case of a symmetric explosion, 
for different values of ^. The probability density for ^ = 10~^ in (a) has been reduced by a 
factor of 100. 

Fig. 5 — Distribution of post-SN systems over dimensionless orbital separations for systems 
that (a) would remain bound (/3 = 0.6), or (b) be disrupted (/3 = 0.4), for different values 

ofC- 

Fig. 6 — Distribution of post-SN binaries over (a) eccentricities, for /3 = 0.6 and ^ = 1.0, and 
over (b) dimensionless orbital separations, for /3 = 0.6 and ^ = 1.0, for different values of c. 

Fig. 7 — Distribution of post-SN binaries over dimensionless circularized orbital separations, 
ac, for different values of /3 and ^. 

Fig. 8 — Distribution of post-SN systems over dimensionless orbital separations of 
circularized orbits, ac, and systemic velocities, Vgyg, for (3 = 0.6 and ^ = 1. 

Fig. 9 — Distribution of post-SN systems over dimensionless systemic velocities, Vsys, for 
Ml = SMq, M2 = 1 M0 (/3 = 0.6), and different values of ^. 

Fig. 10 — Fraction of systems that survive the supernova event, at which kick velocities with 
< Vi >i/2= lOOkms-i arc imparted to the neutron star, as a function of pre-SN orbital 
separation for a 1 M0 secondary and two different masses of the exploding star. In the 
case of a symmetric explosion, systems with Mi = 3.8 Mq remain marginally bound, while 
systems with Mi = 8.6 do not survive. 
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Fig. 11 — Fraction of systems that survive the supernova event as a function of pre-SN 
orbital separation with {dotted line) and without {solid line) an impulse velocity imparted 
to the secondary. 
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